########################################################################
#      File Name:	Figure 5 and S1.R				                             #
#      Date:   	June 2025     			            				               #
#      Author: 	Juan Acevedo-Ossa		                        	         #
#      Purpose:   Replication file for "The Indirect Effects of	       #
#				          Structural Power"                                    #
#      Software: RStudio 2024.09.1                                     #
#      Input Files:  GVC.csv  panel_unga.csv                      	   #	
#      Output File: "Figure 5"  "Figure S1"                            #
#      Machine:  16.0 GB RAM - Intel Core i7 @ 2.50 GHz 16 cores       #
#      Total Processing time:  3 hrs 15 minutes                        #
########################################################################

#1. Install Packages

#install.packages("foreign")
#install.packages("ggplot2")
#install.packages("MASS")   
#install.packages("dplyr")
#install.packages("texreg")
#install.packages("sf")        
#install.packages("spdep")      
#install.packages("spatialreg") 
#install.packages("gridExtra") 

rm(list = ls())


#---------------------------#
# Load required packages
#---------------------------#



library(gridExtra)
library(foreign)
library(ggplot2)
library(MASS)   
library(dplyr)
library(texreg)
library(sf)        
library(spdep)      
library(spatialreg) 


# Get the directory where the script is located
path <- "C:/Users/juanp/OneDrive - City University of New York/The Indirect Effects of Structural Power"
setwd(path)

unga_panel <- read.csv("./Data/panel_unga.csv")
wp <- as.matrix(read.csv("./Data/GVC matrices/GVC.csv"))
np <- nrow(unga_panel) # observations in panel


# Generate unique panel-id name (countrycodes_year) to use in neighbors and weights
unga_panel$id <- paste(unga_panel$countrycodes, unga_panel$year, sep = "_")


# Create a neighbor listw object for pn data
# Note - "nbs_lw_pn" = neighbors_listw_panel
nbs_lw_pn <- mat2listw(x = wp,
                       row.names = unga_panel$id,
                       style     = "W")

# Update w to row-standardized [relevant for marginal effect calculations!]
wp <- listw2mat(listw = nbs_lw_pn)
colnames(wp) <- rownames(wp)

# For LRSS estimates:
n  <- 152                 # Cross-section units
# balanced panel - we have 152 observations per year
w  <- wp[1:n, 1:n]  

# ----------------------------------- #
##  Spatiotemporal lag (STAR) model  
# ----------------------------------- #

unga.lagsartslm <- lagsarlm(formula = alignusimp  ~ lalignusimp + sanct_stock  + ltradegdpuslog  + laidusgdplog  + polyarchy  + pop_log  + loggdppercapita + corruption_control  + imf  + atop + right_wing + dip_visits_usa + unsc_non_p + developing,
                            data    = unga_panel,
                            listw   = nbs_lw_pn)
summary(unga.lagsartslm)

screenreg(unga.lagsartslm)



# ----------------------------------- #
# Import functions
# ----------------------------------- #

source("https://docs.google.com/uc?id=1xmQv4V8BOnXmy0zXpHtulxMbTR4ARx0J&export=download")

#Set a working directory large enough to contain all the simulations

setwd("D:/")


#Simulations take 3.2 hours using a Windows PC 16.0 GB RAM - Intel Core i7 @ 2.50 GHz 16 cores  


mr <- star_mr(model      = unga.lagsartslm,
              lagdv_name = "lalignusimp",
              x_var      = "polyarchy",
              cs_id      = "countrycodes",
              pn_id      = "id",
              sims       = 1000,
              progress   = TRUE)




# ----------------------------------- #
# FIGURE 5
# ----------------------------------- #


ctrys = c("PER", "VEN") 

plots <- list()



for (ctry in ctrys) {
  
  res <- extract_mr(star_sims = mr,
                    shock     = "CHN_2001",
                    response  = ctry,
                    cred_int  = 0.90)
  
  # Cumulative response plot
  plots[[ctry]] <-gplot(data  = res,
                        xlab  = "Time",
                        ylab  = "Response",
                        cumulative = TRUE,
                        title =  paste("Cum resp in ", ctry),
                        color = "blue")
  
  
}


figure5 <- grid.arrange(grobs = plots[1:2], ncol = 2)

figure5

setwd(path)

ggsave("Images/Figure 5.jpg",
       plot = figure5,
       dpi = 300,
       width = 6.26,
       height = 2.56,
       units = "in",
       device = "jpeg")




# ----------------------------------- #
# FIGURE S1
# ----------------------------------- #

ctrys = c("PRY","GTM", "PAN","SLV") 

plots <- list()



for (ctry in ctrys) {

res <- extract_mr(star_sims = mr,
                  shock     = "CHN_2001",
                  response  = ctry,
                  cred_int  = 0.90)

# Cumulative response plot
plots[[ctry]] <-gplot(data  = res,
                      xlab  = "Time",
                      ylab  = "Response",
                      cumulative = TRUE,
                      title =  paste("Cum resp in ", ctry),
                      color = "blue")


}


figures1 <- grid.arrange(grobs = plots[1:4], ncol = 2)

figures1
ggsave("Images/Figure S1.jpg",
       plot = figures1,
       dpi = 300,
       width = 5.07,
       height = 3.85,
       units = "in",
       device = "jpeg")











